Geophysical inversion attempts to find a model of subsurface properties that optimally explains observed data and satisfies geological and geophysical constraints. There are a large number of well known methods of geophysical inversion. These well known methods fall into one of two categories, iterative inversion and non-iterative inversion. The following are definitions of what is commonly meant by each of the two categories:
Non-iterative inversion—inversion that is accomplished by assuming some simple background model and updating the model based on the input data. This method does not use the updated model as input to another step of inversion. For the case of seismic data these methods are commonly referred to as imaging, migration, diffraction tomography or Born inversion.
Iterative inversion—inversion involving repetitious improvement of the subsurface properties model such that a model is found that satisfactorily explains the observed data. If the inversion converges, then the final model will better explain the observed data and will more closely approximate the actual subsurface properties. Iterative inversion usually produces a more accurate model than non-iterative inversion, but is much more expensive to compute.
The most common iterative inversion method employed in geophysics is cost function optimization. Cost function optimization involves iterative minimization or maximization of the value, with respect to the model M, of a cost function S(M) which is a measure of the misfit between the calculated and observed data (this is also sometimes referred to as the objective function), where the calculated data are simulated with a computer using the current geophysical properties model and the physics governing propagation of the source signal in a medium represented by a given geophysical properties model. The simulation computations may be done by any of several numerical methods including but not limited to finite difference, finite element or ray tracing. The simulation computations can be performed in either the frequency or time domain.
Cost function optimization methods are either local or global. Global methods simply involve computing the cost function S(M) for a population of models {M1, M2, M3, . . . } and selecting a set of one or more models from that population that approximately minimize S(M). If further improvement is desired this new selected set of models can then be used as a basis to generate a new population of models that can be again tested relative to the cost function S(M). For global methods each model in the test population can be considered to be an iteration, or at a higher level each set of populations tested can be considered an iteration. Well known global inversion methods include Monte Carlo, simulated annealing, genetic and evolution algorithms.
Unfortunately global optimization methods typically converge extremely slowly and therefore most geophysical inversions are based on local cost function optimization. Algorithm 1 summarizes local cost function optimization.    1. selecting a starting model,    2. computing the gradient of the cost function S(M) with respect to the parameters that describe the model,    3. searching for an updated model that is a perturbation of the starting model in the negative gradient direction that better explains the observed data.Algorithm 1—Algorithm for Performing Local Cost Function Optimization
This procedure is iterated by using the new updated model as the starting model for another gradient search. The process continues until an updated model is found that satisfactorily explains the observed data. Commonly used local cost function inversion methods include gradient search, conjugate gradients and Newton's method.
A very common cost function is the sum of the squared differences (L2 norm) of real and simulated seismic traces. For such a case, the gradient is calculated through a cross-correlation of two wavefields, as shown for the typical full wavefield inversion workflow in FIG. 1. Starting with an estimate of the source wavelet (101), and an initial subsurface model (102), we generate simulated seismic data (103) by propagating waves forward (104) from the source to the receiver locations. The data residuals (105) are formed by subtracting (110) the simulated data from the real seismic data (106). These residuals are then propagated backwards to the subsurface model (107) and cross-correlated with the source wavefield, generated by forward propagation (108) from the source location to the subsurface. The result of this cross-correlation is the gradient (109), on the basis of which the subsurface model is updated. The process is repeated with the new updated model, until the difference between simulated and real seismic data becomes acceptable.
For different cost functions the calculation of the gradient can be different. Still the basic elements of the workflow in FIG. 1 are quite general. The key ideas for the present invention can be trivially modified for cases where alternative cost functions and gradient computations are used.
Iterative inversion is generally preferred over non-iterative inversion, because it yields more accurate subsurface parameter models. Unfortunately, iterative inversion is so computationally expensive that it is impractical to apply it to many problems of interest. This high computational expense is the result of the fact that all inversion techniques require many compute intensive simulations. The compute time of any individual simulation is proportional to the number of sources to be inverted, and typically there are large numbers of sources in geophysical data, where the term source as used in the preceding refers to an activation location of a source apparatus. The problem is exacerbated for iterative inversion, because the number of simulations that must be computed is proportional to the number of iterations in the inversion, and the number of iterations required is typically on the order of hundreds to thousands.
Reducing the computational cost of full wavefield inversion is a key requirement for making the method practical for field-scale 3D applications, particularly when high-resolution is required (e.g. for reservoir characterization). A large number of proposed methods rely on the idea of simultaneously simulated sources, either encoded (e.g. Krebs et al., 2009; Ben-Hadj-Ali et al., 2009; Moghaddam and Herrmann, 2010) or coherently summed (e.g. Berkhout, 1992; Zhang et al., 2005, Van Riel and Hendrik, 2005). Inversion methods based on encoded simultaneous simulation often suffer from cross-talk noise contaminating the inversion result and are commonly limited by the data acquisition configuration (recording data with stationary receivers is a requirement for several of the methods). Methods based on coherent summation typically lead to loss of information. Nevertheless both types of approaches can be very helpful and are the subject of ongoing research.
A different way for reducing the computational cost of full wavefield inversion is by reducing the number of iterations required for convergence, and this is the objective of this invention. The method does not suffer from the typical limitations of the methods mentioned above, but it does not preclude their usage. In fact, it can, in principle, be used in combination with any of the simultaneous-source methods mentioned above, to potentially provide increased computational savings.